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ABSTRACT 


The evolutionary processes generating population genetic structures may differ depend- 
ing on the spatial scale. We examined the population structure of the freshwater snail Radix 
balthica on a local scale in Saanen Valley, Switzerland. We used a combination of a direct 
approach by simulating dispersal with parameters gained out of experiments in an artificial 
pond and indirect molecular genetic approaches by using mitochondrial (COI) and nuclear 
(microsatellite) marker at all sites harbouring R. balthica in this restricted area. 

Contrary to the expectations from the direct methods, populations in homogeneous water bod- 
ies were found to be genetically heterogeneous, independently of geographic distance between 
sites. Instead, several selfing lineages coexisted with little local gene flow among them. 

The population structure observed in the valley could thus not be explained by active dis- 
persal alone. It was rather characterized by independent, probably bird-mediated, colonisation 
events from several sources, in combination with high metapopulation dynamics and a mixed 
mating system preserving this structure. This population structure has rarely been reported 
before, but might nevertheless be typical for passively dispersed, patchily distributed taxa 


(e.g., freshwater invertebrates). 
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INTRODUCTION 


Disentangling the processes governing spa- 
tial genetic structure is essential to understand 
the ecological and evolutionary dynamics of a 
species. Spatial genetic structure is shaped by 
the interaction of demography, natural selec- 
tion, mating system and dispersal-mediated 
gene flow. However, the relative importance 
of these factors may vary according to spatial 
scale. On a local scale, active dispersal may 
play a more important role than long-range 
dispersal and vice versa on a broader scale. 

Furthermore, a mixed mating system can 
substantially affect the distribution of ge- 
netic variation within and among populations 
(Charlesworth, 2003). Predominately selfing 
(sub-) populations are strongly influenced by 
drift due to a very low effective population size 
and are ultimately subject to genetic impover- 
ishment (Jarne & Stadler, 1995). 


“Corresponding author: thaun@bio.uni-frankfurt.de 


219 


A recent study (Pfenninger et al., 2011) on 
the range-wide population structure of Radix 
balthica (Lymnaeidae) revealed a population 
structure characterised by the absence of 
isolation-by-distance, together with rather iso- 
lated and genetically depauperate populations 
compared to the variation present in the entire 
range. This seemingly paradoxical pattern was 
explained by distance-independent, passive 
long-range colonisation and high population 
turn-over rates, low gene-flow and/or strong drift, 
enhanced by the mixed mating system of the 
species (Pfenninger et al., 2011). This pattern 
of population structure, probably typical for pas- 
sively dispersed, patchily distributed taxa (e.g., 
freshwater invertebrates), is likely to impede 
local adaptation (Kawecki & Ebert, 2004). 

Radix balthica occurs in permanent slow- 
flowing rivers and stagnant water bodies, without 
requiring a particular substratum or high water 
quality (Okland, 1990; Glóer & Meier-Brook, 
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1998). The individual active dispersing capacity 
is restricted to the water body the snail lives in. 
The egg clutches are attached to structures, and 
there is no free larval stage allowing dispersal 
(Glöer 8 Meier-Brook, 1998). Still, R. balthica 
is widely distributed over northwestern Europe 
(Pfenninger et al., 2006), implying extensive pas- 
sive transport, as suggested for many freshwater 
taxa (e.g., Frisch et al., 2007). Many gastropods 
have a mixed mating system with outcrossing 
and selfing both common (Jarne & Stadler, 1995; 
Meunier et al., 2001), which is also true for R. 
balthica (Pfenninger et al., 2011). 

To understand the regional and local population 
structure of R. balthica in relation to the range- 
wide structure, we used direct (dispersal studies) 
and indirect (molecular genetic) methods at all 
sites harbouring R. balthica in a restricted area. In 
particular, we asked the following questions: 


What are the expectations for local (intra wa- 
ter body) population structure of R. balthica 
based on active dispersal, and how do they 
correspond to the structure as inferred from 
molecular markers? 

Are the populations connected by gene-flow ona 
regional scale, that is, in an area of 56 km2? 
How does the inferred ‘structure relate to the 
large-scale population over the whole area 

of distribution structure? 


METHODS 
Direct Estimation of Population Structure 


We conducted a series of studies that param- 
eterised computer simulations to estimate the 
life time dispersal based on daily movements. 
This required knowledge of (i) the frequency 
distribution of daily movement distances, (il) 
position preference in the habitat, (iii) depen- 
dence of movement directions on previous 
movement and (iv) dependence of positions 
between individuals (Pielou, 1969). Addition- 
ally, the seasonal water temperature is an 
important factor for the activity of freshwater 
snails and was taken into account. 


Pond Set-up and Recording of Daily Net Move- 
ment Distances: To obtain a distribution of daily 
movement distances, we set up an artificial pond 
of 4 x 4 m in size. The study was conducted 
in June 2010, during which the mean daily 
maximum temperature was 18.8°C. We used 
tap water to fill the pond initially to 20 cm depth 
and inoculated it with 50 litres of natural pond 
water. As substratum, we inserted a 2 cm layer 
of sieved and dried local loam. We let the pond 
settle for two weeks, during which it was naturally 
filled up by rainwater to a depth of 30—35 cm and 
accumulated sufficient debris as a food source 
for the snails. Thirty individually marked adult 
snails > 10 mm from a laboratory stock popula- 
tion were released in the middle of the pond. 
After an acclimatisation period of one week, we 
recorded the position of living individuals every 
24 hours relative to one corner of the pond for 
a period of 21 days. The daily net dispersal 
distances between two records were calculated 
as straight lines. We also recorded the dispersal 
directions and their changes between successive 
days. Individuals not surviving the measurement 
period were discarded from the analysis. The dis- 
tribution of daily net movements over the study 
period was fitted to an exponential function. 


Testing for a Preferred Position in the Habitat: 
Observations in the field and during the study 
suggested that the snails tended to stay close 
to the rim. We therefore measured the distance 
of each snail to the closest shoreline on the last 
day of the survey, to compare the observed 
distribution with that expected under a com- 
puter simulated movement model. For this 
simulation, a virtual snail was set at a random 
position within a virtual 4 x 4 m pond. A random 
variate was drawn from the observed daily 
net movement distribution and assigned to a 
random direction. If a snail would have crossed 
the rim of the pond with a movement in the 
simulation, it was reflected at the rim with the 
opposite direction and moved on the remaining 
distance. This was repeated 21 times for each 
snail and after the last move, the distance to 
the closest rim was recorded. The simulation 
was run 1,000 times for 29 snails. 


—> 


FIG. 1. Spatial distribution of the sampling sites along the Saane valley in Switzerland (square) and 
three known neighbouring populations (red dots). The position of the 18 sampling sites within the water 
bodies (labelled A-D) are indicated by coloured circles representing the fraction of the mtDNA lineages 
(a, B, y), together with the cyto-nuclear composition of the individuals sampled. Vertical bar represents 
a single individual. Column shows the estimated nuclear composition of the multilocus microsatellite 
genotype (dark green = cluster |, light green = cluster II, yellow = cluster III, blue = cluster IV) as indi- 


cated by STRUCTURE. 
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Test for Independence of Previous Moves 
and of Individual Positions: Independence of 
movement directions from previous moves 
was tested by comparing the number of 
changes in consecutive movement direc- 
tions to random expectations with a x° test. 

To test whether the snail's positions (and thus 
the movements) were independent from each 
other or whether they tended to aggregate or 
avoid each other, we recorded the average 
distance between the snails after 21 days 
and compared it to the expected distribution 
under a movement simulation model. We 
used basically the same simulation model as 
described above with two modifications. First, 
as previous results indicated that the snails 
moved almost exclusively along the rim of the 
pond, we considered only the distance vec- 
tors parallel to the rim from the observed daily 
net movement distribution. Second, because 
consecutive moves were not independent 
from previous moves, direction changes were 
randomly assigned according to the observed 
direction-change probability distribution. The 
average distance along the rim between 
all snails was recorded after 21 days. The 
simulation was run 10,000 times for 29 snails. 


Temperature Dependence of Activity: To infer 
the temperature dependence of activity in Ra- 
dix, individually marked snails (n = 10) of each 
of two laboratory stock populations were put 
in a 30 litre glass aquarium filled with aerated 
water in a climate chamber (one aquarium for 
each stock population). Food ad libitum (pow- 
dered TetraMin fish-food, Tetra) was provided. 
We increased the ambient temperature of the 
chamber from 2°C to 34°C in steps of 2°C/h. 
When the next temperature level was reached, 
which took about one hour, the temperature 
was held constant for the measurement of 
individual movements for one hour. Total dis- 
tances covered in this hour by each individual 
were recorded. Movement per hour was aver- 
aged over both investigations. To establish 
the relationship between activity and water 
temperature, we calculated a linear regression 
between movement (dependent variable) and 
temperature (independent variable). 


Modelling Life-Time Dispersal. The results of 
the simulations described above were used 
to design and parameterise a life-time disper- 
sal model for Radix. lt was assumed that the 
snails have a life span of approximately 260 
days. This matches field observations of an 


annual life cycle (Dillon, 2000) and equals the 
number of the days in the River Saanen Valley 
with a mean temperature above 5°C, exclud- 
ing the time for hibernation. To simulate the 
life time dispersal distance of a snail in a two 
dimensional habitat of unrestricted size, it was 
initially placed at position zero. A random variate 
from the exponential function of the observed 
daily movement distance distribution along the 
shore line was drawn. To take the temperature 
dependence of movement activity in Radix into 
account, the assumed life time was divided into 
eight periods of 30 days (corresponding to the 
months March to November) and two shorter 
periods of ten days at the beginning and the 
end (each during February and December). 
The determined temperature-dependent activity 
function was then used to determine the correct 
drawn movement distance according to the av- 
erage monthly mean temperature in the Saanen 
Valley. Climate data for the Saanen Valley were 
taken from http://www.meteoschweiz.admin. 
ch. Whether the direction of movement for the 
subsequent day changed or not was randomly 
determined according to the observed direction 
change probability. This procedure was repeat- 
ed 260 times and the distance of the final posi- 
tion to the starting point of the snail recorded. 
We simulated 100,000 snails in total. 

Wright (1946) calculated the neighbourhood 
size for a linear-distributed population in a two- 
dimensional habitat as N = V2 p - ø - d , with o 
being the variance of dispersal along the habitat 
from parents relative to their offspring at the 
same developmental stage and d as the ob- 
served density of breeding individuals. The strip 
length where parents potentially can originate 
from and no population differentiation can be ex- 
pected thus equals about 3.55 o (Wright, 1946). 


Indirect Estimation of Population Structure 


Sampling Design: Sampling was conducted in 
a part of the river Saanen Valley, Switzerland. 
We searched for Radix populations in all natural 
and artificial water bodies in the valley between 
Montbovon (46.628°N, 7.044°E) and the Lac 
de la Gruyère (46.623°N, 7.101°E) and below 
1,000 m altitude. Two hundred and seventy 
individuals from 18 sites were sampled. Radix 
were found in only four out of 35 water bodies: 
two small creeks, a gravel pit and a reservoir 
of the Saanen river (Table 1, Fig. 1). Along the 
creeks where Radix occurred, all adult snails 
(> 8 mm) that could be found were sampled 
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along a stretch of 5 m, every 50 m. In the gravel 
pit, all adult snails were sampled and pooled 
per 10 m shoreline. If the density was high, a 
maximum of 22 adult snails were taken per 
section. At the reservoir, all accessible parts 
of the shoreline were searched and 11—30 
snails were sampled at places where Radix 
were present. Individuals were immediately 
preserved in 100% ethanol. The distances 
between sampling sites were measured along 
the shoreline (Table 2) for comparison to the 
neighbourhood size calculated from the simula- 
tion of lifetime dispersal. 


Microsatellite and Mitochondrial Haplotype 
Analysis: DNA was extracted using a glass 
fibre DNA extraction method following the pro- 
tocol of Ivanova et al. (2006). About 1 mm’ of 
alcohol-preserved snail tissue was mixed with 
vertebrate lysis buffer containing 100 mM NaCl, 
50 mM Tris-HCI (pH 8.0), 10 mM EDTA (pH 
8.0), 0.5% SDS and 20 mg per ml proteinase 
K. Samples were incubated over night at 56°C. 
After digestion, samples were mixed with bind- 
ing buffer (BB), composed of 6 M guanidinium 
thiocyanate, 20 mM EDTA, 10 mM Tris-HCl 
(pH 6.4), 4% Triton X-100, transferred into a 
Pall1 glass fibre plate (GF plate, AcroPrep™), 
placed on top of a 2 ml square well block and 
centrifuged to bind DNA to the GF membrane. 
DNA washing was carried out in two steps: (1) 
Protein wash buffer: BB and EtOH 96%, (2) 
Wash buffer: 60% EtOH, 50 mM NaCl, 10 mM 
Tris- HCI (pH 7.4), 0.5 mM EDTA (pH 8.0). After 
washing, residual ethanol was evaporated by 
incubating the GF plate at 56°C for about 30 
minutes. DNA was eluted by incubating each 
GF membrane with 50 ul of prewarmed (56°C) 
ddH,0 followed by centrifugation using a 96- 
well collection plate. 

Cytochrome oxidase subunit | (COI) frag- 
ments were amplified using polymerase chain 
reaction (PCR) performed with Invitrogen 
Taq DNA polymerase and universal primers 
(Folmer et al., 1994). Sequencing reactions 
were performed using the ABI Prism Big Dye 
terminator kit (Perkin-Elmer, Norwalk, CA, 
USA). Sequenced fragments were read on an 
ABI Prism 3730 automated capillary sequencer 
(Applied Biosystems). 

All snails were genotyped for eight highly 
polymorphic microsatellite loci (Salinger & 
Pfenninger, 2009). Multiplex Microsatellite 
amplification was carried out using the Quiagen 
Type-it™ microsatellite PCR Kit with fluores- 
cent dye labeled forward primers (Salinger 
& Pfenninger, 2009). PCR products were 


separated and visualized on an ABI Prism 
3730 sequencer (Applied Biosystems) with 
GeneScan™ 500-LIZ™ size standard (Applied 
Biosystems). Microsatellite allele lengths were 
recorded using GENEMAPPER 4.0 software (Ap- 
plied Biosystems). 


Genetic Variability: Mitochondrial sequences 
of the COl-fragment of 482 bp length were 
aligned using ClustalW (Thompson et al., 
1994) as implemented in MEGA 4.0 (Tamura 
et al., 2007) and corrected by eye. We used 
this alignment to calculate a parsimony hap- 
lotype network in TCS 1.21 (Clement et al., 
2000). We nested haplotypes in a hierarchical 
fashion as described by Pfenninger & Posada 
(2002). Identification of the sampled snails as 
R. balthica was done by comparing the respec- 
tive barcodes to the reference sequences of 
Pfenninger et al. (2006; GenBank Accession 
no. DQ980030—DQ980193) using a neighbour- 
joining tree and monophyly with the previously 
identified haplotypes. This was necessary, 
because other Radix species are present in 
the region and specific identification is only 
possible with molecular markers (Pfenninger 
et al., 2006). Mitochondrial genetic distances 
among the water bodies were calculated as 
Neis standard genetic distances (Dr; Nei, 
1972) with MEGA 4.0 (Tamura et al., 2007). 
We used GenePop (Raymond & Rousset, 
1995a) to determine allele and genotype 
frequencies of the microsatellite markers, for 
the different sampling sites, the four water 
bodies and the entire sample. To investigate 
population differentiation, Fs; was calculated 
using FSTAT 2.9.3.2 (Goudet, 1995). FSTAT 
was also used to estimate the expected (He) 
and observed heterozygosities (Ho) for each 
sampling site and to test for deviation from 
Hardy-Weinberg Equilibrium among sampling 
sites within water bodies (1,000 permutations, 
two-sided test). Bonferroni correction for mul- 
tiple tests was applied to each of these tests. 
Nuclear genetic distances among the water 
bodies (Dc; Nei, 1972) were calculated using 
Genetix 4.04 (Belkhir et al., 1996—2004). The 
microsatellite data were also subjected to a 
factorial correspondence analysis with Genetix 
4.04 (Belkhir et al., 1996—2004) to visually ex- 
plore the distribution of genetic variation. 


Sensitivity Test: Because our data set contained 
sampling sites with relatively few individuals, 
we conducted two initial simulation studies 
to assure that observed deviations from the 
null hypothesis of no population structure are 
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not due to sampling artefacts. First, we simu- 
lated 300 times a truly panmictic population of 
200,000 individuals with a diploid genotype of 
eight independent loci with the same possible 
number of allelic states as observed in our 
sample of natural populations, respectively, 
using EASYPOP (Balloux, 2001). The muta- 
tion rate was set to 0.05 and a model with 
mutations egually likely to occur at any of K 
possible sites (KAM) was chosen. After 1,000 
generations, a dataset from the final genera- 
tion of 235 individuals was generated, divided 
into 18 groups ofthe same size as our sample 
of natural populations. A second dataset was 
generated almost identically according to the 
sampling parameters at sampling site A. The 
number of individuals was set to 77 and the 
number of groups was set to 7. 

The resulting datasets were tested for devia- 
tion from the null hypothesis of no population 
structure with an AMOVA (Analysis of Molecular 
Variance, 10,000 replicates) in Arleguin 3.11 
(Excoffier et al., 2005). 


Inference of Population Structure: Assignment 
of individuals to genotype clusters without 
a priori classification was performed with 
STRUCTURE (Pritchard et al., 2000). All MCMC 
runs were repeated five times for each value 
of K (K = 1 to K = 20) for 1,000,000 genera- 
tions with 50,000 burn-in steps, applying the 
admixture model. We also used a new model 
developed for Structure that integrates in- 
formation about sampling locations. This can 
be useful if the cluster identification without a 
priori information would not suffice to detect 
clusters in a dataset (Hubisz et al., 2009). We 
used the same settings as in the model without 
the location-prior. Both models were analysed 
using the maximum LnP(D) value. LnP(D) is 
the log likelihood of the observed genotype 
distribution in K clusters. 

Fs; among water bodies was calculated with 
Genetix 4.04 (Belkhir et al., 1996-2004). An 
exact test of population differentiation (Ray- 
mond & Rousset, 1995b) was performed using 
Arlequin 3.11 between all pairs of sampling 
sites within the respective water bodies, with 
1,000,000 steps of Markov Chain length and 
1,000,000, dememorization steps. 


Isolation by Distance: We tested for Iso- 
lation-by-Distance among sampling sites 
using a Mantel’s test (Mantel, 1967). A cor- 
relation between Fsr / (1 — Fsr) and the 
natural logarithm of pair-wise geographical 
distances was obtained using the ISOLDE 


program of the GENEPOP package (Ray- 
mond & Rousset, 1995a; Rousset, 2008). 


Estimation of Selfing Rates: Genotype frequen- 
cies at codominant marker loci may convey 
information on mating systems (David et al., 
2007). The classical approach is to measure 
heterozygote deficiencies and obtain the selfing 
rate s from Fıs = s / (2-s), assuming inbreeding 
equilibrium. A drawback of this estimate is that 
heterozygote deficiency may not only be due 
to selfing. Biparental inbreeding and partial 
dominance as well as such technical artefacts 
as null alleles also increase the proportion of 
observed homozygotes. Therefore, in addition 
to this classical approach, we used a method 
implemented in RMES (David et al., 2007) to 
derive estimates of s(g>), based on a point 
estimate of the second-order heterozygosity 
disequilibrium and a maximum likelihood of 
the whole distribution based approach of Smu). 
These estimates are independent of Fıs and less 
influenced by technical biases. 


Assignment Tests: To compute the probability 
of each individual being derived from external 
populations or being a resident (i.e., not a first- 
generation migrant) in the Saanen Valley, we 
used GENECLASS2 (Piry et al., 2004). As a ref- 
erence dataset we used 80 populations, covering 
the approximate species range, from a previous 
study on R. balthica (Pfenninger et al., 2011). 
This software computes genetic assignment 
statistics to assign or exclude reference popula- 
tions as the origin of individuals on the basis of 
multilocus genotype data. It is thus also possible 
to obtain direct estimates of dispersal through 
the detection of immigrant individuals (Rannala 
8 Mountain, 1997; Paetkau et al., 2004). For the 
assignment test and the test for first generation 
migrants we selected the likelihood ratio L home 
/L_max and a Bayesian approach by Baudouin 
& Lebrun (2001) with an assignment threshold of 
scores of 0.01 and tested for probabilities using 
a simulation algorithm by Paetkau et al. (2004) 
with 100,000 simulated individuals and an alpha 
error set to 0.05. This approach provides a prob- 
ability estimate of exclusion from or inclusion 
in a population, given only the genotype of the 
animal in question, and estimates of the allele 
frequency of the population in question. This 
simulation approach is sensitive to sample size 
and the degree of genetic differentiation among 
the populations. The significance of results 
was estimated with Monte Carlo re-sampling 
techniques as implemented in GENECLASS2 
(Paetkau et al., 2004). 
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RESULTS 
Movement Dispersal Patterns 


During the 21 days of the study, 402 moves 
were recorded from 27 snails. Daily distances 
ranged from 0.01 to 3.86 m, with a median of 
0.43 m. Distances covered in 21 d ranged from 
3.22 m to 23.02 m, with a mean of 13.60 + 4.96 
m (mean + s.d.). The frequency distribution of 
daily dispersal distances fitted a function with 
exponential decay (y = 85.7 e-0.23x: 72 = 0.88, p 
< 0.001, n = 1,000). 

The observed average distance between indi- 
viduals after 21 days was 5.25 += 3.88 m (mean 
+ s.d.). Average inter-individual distance from a 
random movement simulation was 5.53 + 3.95 
m (mean t s.d.). The observed value falls well 
within the 95% confidence interval. The null 
hypothesis of random distribution could thus 
not be rejected. 

Movement along the rim was preferred over 
movement away from the rim (x2 = 47.12; n 
= 395; df. = 3; p < 0.001; comparison with 
a uniform distribution). The snails lived and 
moved almost exclusively close to the rim of 
the artificial pond. Their observed distribution 
relative to the rim was significantly different 
from the. expected distribution obtained from 
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a random movement simulation. Directions of 
second moves were dependent on direction of 
the first moves (x2 = 22.69; n = 365; d.f. = 2; p< 
0.001; comparison of deviation from first move 
with a uniform distribution). Snails maintained 
directions on consecutive days in 62% ofcases. 
Random movement could thus with respect to 
direction not be assumed. 

Movement of R. balthica was strongly tem- 
perature dependent. A linear regression for 
the distance moved per hour against water 
temperatures was highly significant (Fig. 2). 

Computer simulated mean effective lifetime 
dispersal in a two-dimensional habitat was 
48.25 m (s.d. = 35.75 m) in 260 d. Minimum 
simulated effective lifetime dispersal distance 
was 0.11 m and maximum distance was 
229.62 m. 


Direct Estimation of Population Structure 


The standard deviation of life time dispersal 
was 35.75 m. The spatial extent of a neigh- 
bourhood was thus 126.9 m of shore line. If 
we compare this to the distances between 
the sampling locations (Table 2), we do not 
expect population differentiation at least among 
sampling points in water body A and between 
sampling points D2 and D3. 


20 25 30 35 


Temperature [°C] 


FIG. 2. Linear regression of R. balthica movement per hour against water temperature. Movement per 
hour are averages measured from 10 snails. The regression was highly significant (p < 0.0001). 


TABLE 2. Distances in meters between the sampling points along the shoreline within the water bodies. Asterisks show significance level of an exact test 
of population differentiation (*p < 0.05; **p < 0.01; ***p < 0.001) between all pairs of samples within the water bodies. 
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Incidence of Radix balthica in the Sampling 
Region 


Despite the presence of many seemingly 
suitable water bodies in the Saanen Valley, 
Radix snails were found only in a small fraction 
of them (4 out of 35). These tended to be the 
larger ones, although the species was absent 
from one ofthe larger gravel pits in the immedi- 
ate vicinity of water body A. 


Mitochondrial Haplotype Variability 


At least 482 bp of the COl-gene fragment 
were obtained and manually aligned for 236 
individuals from the Saanen Valley. All were 
identified as R. balthica by comparing the se- 
guences to the reference seguences (data not 
shown). There were 28 polymorphic sites in our 
sample and 19 haplotypes (GenBank, acces- 
sion no. GU735965-GU736200). Haplotypes 
were designated according to the nomenclature 
of Cordellier & Pfenninger (2009). Statistical 
parsimony analyses clustered these haplotypes 
into three lineages (a, B, y, Fig. 3). Lineage a 
included 9 haplotypes from more than half the 
individuals analyzed (n = 169), B comprised 4 
haplotypes from 23 individuals and y included 6 
haplotypes from 34 individuals. The distribution 


of haplotype lineages among sampling sites is 
shown in Figure 1. Genetic diversity within lin- 
eages was between 0.4 and 5.9% and genetic 
divergence among lineages ranged from 0.1 to 
7.1%. The genetic diversity at each sampling 
site and within each water body was measured 
by two estimators, haplotype diversity (Dhap), 
which ranged from 0.08 to 0.8, and genetic 
diversity (Dt), which ranged from 0 to 14% 
between sampling sites with a mean of 5% 
differentiation (Table 1). 


Nuclear Microsatellite Variability 


Two-hundred-and-thirty-five individuals were 
successfully genotyped at 8 microsatellite 
loci. The number of alleles per locus ranged 
from 4 to 16, with an average of 4.4 (Table 1). 
The mean number of alleles (N,) ranged from 
1.85—5.13; Nei’s mean genetic diversity had 
values from 0.2 to 0.59 with an overall mean 
of 0.39. All microsatellite loci showed consis- 
tent heterozygote deficiency with an overall 
observed heterozygosity of 21.4% and an 
expected heterozygosity of 48.6% for the sam- 
pling sites. Values of expected heterozygos- 
ity (He) ranged from 0.19—0.58 and observed 
heterozygosity (H,) from 0.1—0.36. Detailed 
values are shown in Table 1. Significant devia- 


FIG. 3. Minimum spanning network of R. balthica. Each circle represents a haplotype, its size being 
proportional to the frequency of a certain haplotype. Connecting lines represent one mutational step. 
Small black circles depict haplotypes that where not sampled. The lineages are named a, B, v. 
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TABLE 3. The grouping of the nuclear genotypes 
into four clusters (I, II, HI and IV) is based on the 
K from the Srrucrure analysis. Lineage a, B and 
y refer to the inferred mitochondrial lineage. The 
Rxc exact contingency table test on independence 
was significant (d.f. = 6, x2 = 16.6, p < 0.001). 
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tions from Hardy-Weinberg equilibrium were 
observed in nearly half the loci in the sampling 
points, even after Bonferroni's correction. We 
attribute this systematic deviation to a mixed 
mating system common in pulmonates (Jarne 
& Stadler, 1995) and described for R. balthica 
(Salinger & Pfenninger, 2009). 


Cyto-Nuclear Association Test 


There was a weak, albeit significant associa- 
tion of the two nuclear clusters (inferred from 
the AK) with mitochondrial haplotype lineages 
(N = 196, p < 0.001, x2 = 16.6, d.f. = 6, exact 
4x3 Rxc contingency test). Slightly less indi- 
viduals of nuclear cluster | contained lineage 
y haplotypes than expected by chance. The 
same was true for nuclear cluster III individuals 
with lineage a haplotypes (Table 3). 


Sensitivity Test 


The simulated datasets indicated that the 
actual number of loci used, their variability 
and the applied sampling design were capable 
of correctly keeping the null hypothesis of no 


population differentiation if it is true. This was 
true in 291 out of 300 tested cases, for the full 
dataset. For the reduced dataset, equivalent to 
water body A, this was true in 287 out of 300 
tested cases, corresponding in both cases to 
a type | error of less than 0.05. 


Inference of Population Structure 


Applying a location prior to the STRucTURE 
analysis found the highest mean likelihood 
value for K = 4. Seventy individuals were at- 
tributed to the first cluster (1), 51 to the second 
(Il), 22 to the third (III), 53 to the fourth (IV). 
However, even though the number of clusters 
found corresponded to the number of water 
bodies examined, there was no clear associa- 
tion between them (Fig. 1). 

The overall Fst among all sampling sites was 
0.15, ranging from 0.08 to 0.25. Results of the 
test for exact population differentiation within 
the water bodies showed highly significant dif- 
ferentiation between several sampling points 
(Table 2). 

We found no pattern of overall isolation-by- 
distance (r2 = 0.0005, Mantel test p = 0.491). 


Population Assignment Test 


Applying the simulation approach, 213 out of 
235 animals (91%) were genetically assigned 
as residents rather than immigrants. Assign- 
ment was probable to three outside populations 
for 13 individuals to a neighbouring population 
in northern Switzerland (AUG — 46.615°N, 
7.181°E) and seven and two individuals to 
either of two populations in southern Sweden 
(SEM — 57.155°N, 16.419°E; EDA — 59.834°N, 
17.878°E). 

Testing for first generation migrants highlight- 
ed a single individual with most likely external 
origin in AUG. Eight individuals identified as first 


TABLE 4. An analysis of molecular variance (AMOVA) for microsatellite data for the 4 clusters 


estimated with STRUCTURE (K = 4). 


Sums of Variance % of 
d.f. squares components variation p-value 
STRUCTURE Clustering = 4 
Among water bodies 3 79.1 0.23 13.56 0.001 
Within water bodies 466 690.5 1.48 86.44 0.001 
Total 469 769.6 1.71 
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generation migrants, however, were assigned 
to other sampling sites within the study. Four 
individuals from A were assigned to C, both 
one individual from B was assigned to D and C 
respectively and three individuals from D were 
assigned to C. 


Estimation of Selfing Rates 


The inbreeding coefficient Fis of the sampling 
sites ranged from 0.31 to 0.57 with an overall 
mean of 0.51. The selfing rate s estimated from 
Fis ranged from 0.02 to 0.8, the overall mean 
was 0.34. The selfing rate derived with RMES 
(David et al., 2007) differed from the former 
estimate not so much in the mean, ratherin the 
variance, as s(g2) ranged from 0 to 0.9 with an 
overall mean of 0.33. The Maximum Likelihood 
estimate of s(y,) ranged from 0.01 to 0.84 with 
an overall mean of 0.15. 7 


DISCUSSION 


Observed Local Population Structure Does 
Not Match Expectations from Active Dispersal 


Estimating dispersal from observational data 
(Cameron & Williamson, 1977) or in combina- 
tion with computer simulation has been suc- 
cessfully applied before in several studies on 
land snails (Baur & Baur, 1993; Pfenninger et 
al., 1996). These studies showed that simu- 
lated dispersal was similar to dispersal actually 
measured in natural habitats which justified its 
application in our study. 

From active dispersal alone, we would have 
expected a homogeneous population within 
water body A and isolation by distance in water 
body C because: (i) the estimated neighbour- 
hood where no differentiation may be expected 
comprises ~130 m of shoreline, which ap- 
proximately corresponds to the shoreline of 
water body A, and (ii) our marker system would 
probably not be able to distinguish even a mul- 
tiple of this measure. Nevertheless, sensitivity 
analysis of our data proved their suitability for 
correctly supporting the null hypothesis of no 
differentiation. However, the population differ- 
entiation inferred by nuclear and mitochondrial 
markers was much larger than expected (Table 
4). There are several possible explanations for 
these discrepancies. First, the direct estimate 
might be exaggerated due to our study setup 
because (i) natural water bodies with diverse 
habitat structure might reduce the actual 


neighbourhood size range due to preference 
of certain spots (Bovbjerg, 1968), and (ii) the 
density in the artificial pond might have been 
too high, which may have lead to crowding 
effects. However, since snails in most water 
bodies occurred in much higher abundances 
than in the study and were uniformly distributed 
along the shoreline (personal observation MP), 
these two issues should be of minimal impor- 
tance. Second, the parameters from which 
neighbourhood size was calculated could have 
been derived from incorrect assumptions. For 
example, measurement of the temperature 
dependence of activity might not translate in a 
linear fashion to the field. Snails were kept in dif- 
ferent temperature regimes for only one hour, so 
the metabolism might have lagged behind the 
actual temperature or the recorded distances 
might just show acclimatisation activity (Bailey 
& Johnston, 2005); for example, the attempt 
to escape high water temperatures that would 
have seized after a while. However, the linear 
response to temperature shifts argues against 
this latter possibility. Rather the relative correc- 
tion of the dispersal distances measured in the 
field avoided the introduction of a systematic lag 
bias. Also a non-linear dependence of air tem- 
peratures (as reflected in the climate data) to 
water temperature would have biased the simu- 
lation results. However, except for water bodies 
influenced by geothermic or anthropologic activ- 
ity, temperatures of surface water bodies, such 
as those observed here usually correspond in 
a linear fashion to air temperature. 

The use of a Radix balthica population from 
outside the valley for the dispersal study might 
have biased the estimate if dispersal activity is 
population specific (Terblanche et al., 2009). 
Heritable variation in dispersal activity cannot 
be excluded; however, it is highly unlikely that 
differences on the order of magnitudes neces- 
sary to explain the discrepancy exist among 
the populations used for the direct approach 
and those observed in the field. 

Third the mixed mating system could influ- 
ence the population structure. High selfing 
rates, which ranged from 0 to 0.9 (mean 0.33; 
Table 1) in our study, should increase local 
inbreeding, thereby strongly decreasing the 
actual neighbourhood range (especially in 
the case of complete selfing) and thus lead 
to stronger spatial differentiation. However, 
in this case we would expect a pronounced 
isolation-by-distance pattern within water 
bodies at dispersal-drift equilibrium, which is 
not the case. 
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To explain the observed population structure 
based on neighbourhood range alone, the 
neighbourhood must not exceed a few square 
meters, which seems very unlikely, given the 
observational data. Even though we cannot be 
sure that we did not overestimate the neigh- 
bourhood range, lack of active dispersal is 
probably not sufficient to explain the structure 
observed. A non-equilibrium situation therefore 
remains as the most probable explanation. 

Independent colonisation of the valley from 
several sources in combination with metapopu- 
lation dynamics and the mixed mating system 
may be provide the processes causing the 
observed pattern. The haplotype and nuclear 
diversity found in the valley encompassed a 
substantial fraction of the diversity in the entire 
species range (Cordellier & Pfenninger, 2009), 
indicating colonisation of the valley from differ- 
ent sources. Indeed, the deviation of observed 
data from expectations results mainly from the 
presence of different nuclear clusters/mito- 
chondrial lineages at certain sampling points 
within water bodies (Fig. 1). For example, sam- 
pling point A7 is differentiated from the other 
sampling points in water body A. This part is 
inhabited by a different lineage, indicating two 
independent colonisation events. 

The weak albeit significant association of 
mitochondrial lineages with nuclear clusters 
indicated that there is either high ongoing gene- 
flow, constantly supplying individuals from the 
same populations or that the mixed mating 
system retards complete mixing of individuals 
from several immigration events in the past. 
The low estimate of first generation migrants 
both from within the valley and from the entire 
species range argues for the latter scenario. 

Absence of the species from many appar- 
ently suitable water bodies during the sampling 
program, but subsequent colonisation of some 
(personal observation MP and TH) indicates 
a substantial metapopulation dynamic, which 
may influence population structure strongly. 
Metapopulation dynamics with frequent extinc- 
tion and recolonisation events are common 
in other freshwater snails (Viard et al., 1997; 
Charbonnel et al., 2002; Facon & David, 2006). 
Since Radix shells are very thin-gauge, it is not 
possible to find empty shells which could give a 
hint at past occupations from other water bodies 
in the valley. | 

While the snails are not able to cross land bar- 
riers, passive dispersal must be assumed. For 
freshwater snails in particular, waterfowl have 
been suggested as dispersal vectors (Malone, 


1965; Boag, 1986; Bilton et al., 2001). This is 
even more likely since land snails can survive 
passage through a bird's gut (Wada et al., 
2011) and freshwater snails passage through 
a fish's gut (Brown, 2007). However, direct 
data on the gualitative and guantitative impact 
of bird-mediated dispersal is lacking and thus 
urgently needed to fully understand the popula- 
tion dynamics of freshwater invertebrates. 
The fine-scale population structure here in- 
ferred by nuclear and mitochondrial markers 
thus reflects the species range-wide structure 
(Pfenninger et al., 2011), shaped by processes 
such as long-range dispersal, large-0scale 
metapopulation dynamics and low levels of local 
passive dispersal. Active dispersal had, even in 
smallest water bodies, a nominal impact on local 
population structure. Furthermore, the hereby 
formed structure is preserved by the mixed 
mating system. We found no evidence that local 
processes to influence the observed structure. 
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